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The genus Rana Linnaeus, 1758 includes Western Palearc- 
tic brotvn frogs; within this group eight species occur in 
Europe and 32 in Asia (Speybroeck et al. 2016, Yuan et 
al. 2016). The group underwent a basal post-Messinian ra¬ 
diation about 4 million years ago when they split into five 
major clades (Veith et al. 2003). Climatic changes with 
temperature oscillations that followed the Messinian cri¬ 
sis, restricted most brown frog species in their refugia with 
enough time to speciate, although the diversification into 
major clades was completed before the first Pleistocene 
glaciations (Veith et al. 2003). One of the three mono- 
typic lineages that arose during the diversification, is the 
Stream Frog - Rana graeca Boulenger, 1891. Picariello 
et al. (2002) used Si satellite DNA as a taxonomic marker 
to raise the former subspecies R. graeca italica (occurring 
in the Apennines) to the species level (Italian Stream Frog 
- R. italica Dubois, 1987), which restricted the distribution 
of R. graeca to the Balkans and Peloponnese in a continu¬ 
ous form (Fig. lA). Phylogenetic analysis conducted so far 
(Veith et al. 2003, Yuan et al. 2016) show that the clos¬ 
est relatives to R. graeca are R. latastei Boulenger, 1879 
and R. dalmatina Fitzinger, 1839 whereas the, apparently 
closest, R. italica revealed to be even a distinct evolution¬ 
ary lineage (Veith et al. 2003). Afore mentioned Rana spe¬ 
cies are all Mediterranean endemics with the exception of 
R. dalmatina that is widely distributed over Europe (Spey¬ 
broeck et al. 2016). The distribution extent of R. graeca 
might have been wider than today before the glaciations, 
since fossils are found outside their current distribution 


(Veith et al. 2003). Nowadays, in regard to other Rana 
species, R. graeca occurs in sympatry only with R. dalma¬ 
tina and R. temporaria Linnaeus, 1758. The most mark¬ 
ing characteristic that discriminates it from the latter two 
species is a mottled throat that has a pale central area re¬ 
sembling a white line dividing the marbled throat pattern 
in two equal parts (Lelo 2013). Rana graeca is specialized 
to canyons of fast flowing rivers, mountain streams and 
brooks (Speybroeck et al. 2016). More rarely, the species 
can be found near closed mountain lakes and estuaries or 
river deltas not far from the seaside (Sunje et al. 2018). The 
vertical distribution of the Stream Frog goes from 200 to 
2,000 m a.s.l. (Dzukic 1968). 

Although the phylogenetic position of R. graeca is well 
studied in respect to other Rana species (Veith et al. 
2003, Yuan et al. 2016), evolutionary relationships within 
R. graeca populations are completely unknown. This work 
aims to provide new insights to answer this question by 
analyzing a partial sequence of the cytochrome-b (cyt-b) 
of samples collected along the species range. 

Sampling was done in the field and from the material de¬ 
posited in various collections (Natural History Museum of 
Vienna, Hungarian National History Museum in Budapest, 
and private Zoological collection “Lelo” in Sarajevo). Field 
work was conducted in Bosnia and Herzegovina (BIH) and 
Montenegro (Table 1). Small toe clips were sampled and 
preserved in 95% ethanol. The material from the Zoologi¬ 
cal collection “Lelo” was stored in formalin, whereas the 
museum material was stored in ethanol prior to taking tis- 
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sue samples. From these collection samples, a piece of tis¬ 
sue from the hind leg was taken and preserved in 95% etha¬ 
nol. In total, 22 samples, well spread over the entire range of 
R. gracea, were included in this study (Fig. iB, Table 1). Ad¬ 
ditionally we obtained a sample of a R. dalmatina from the 
location Vrelo Bosne (VB, Table 1) during field work in BIH. 

Total genomic DNA was extracted manually from 
the toe clips using the Qiagen DNeasy Blood and Tissue 
Kit (Quiagen) and from the museum samples using the 
QIAamp DNA Micro Kit (Qiagen), following the manu¬ 
facturer’s instructions. The samples taken from museum 
specimens were kept in 1 ml of PBS buffer for 24 hours pri¬ 
or to DNA extraction to enhance the purity of tissues. The 
quality and amount of extracted DNA was measured using 
a NanoDrop machine (Thermo Fisher Scientific) and after¬ 
wards diluted to a working concentration of long/pi. 

Partial sequences of cyt-b were amplified using primer 
sequences (Rana-Cytb-F2, and Rana-Cytb-R2) and tem¬ 
perature profiles given in Vences at al. (2013). The cyt-b 
was amplified in a total volume of 20 pi which included 
0.8 pi of 0.2 pM of each primer, 2 pi of each of the: 200 pM 


of dNTP mix, 1 x Buffer, and 2.5 mM MgCL Additionally, 
0.08 pi of 0.5 U of Taq DNA polymerase, and 7.3 pi of MiliQ 
Water were finally added to 5 pi template DNA. For collec¬ 
tion samples, 2 pi of 1 x Buffer containing MgCl^ (instead 
of separately adding MgCy was used for the PCR reaction 
together with 2 pi of 1 x BSA buffer. For all amplifications, 
contamination was checked using negative controls on gel 
electrophoresis. 

Before sequencing, PCR products were purified using 
the high pure PCR product purification kit (Roche) follow¬ 
ing the manufacturer’s protocol. DNA sequences of forward 
strands were obtained for all PCR products following the 
ABI Prism Big-Dye Terminator Cycle Sequencing standard 
protocol. The sequencing reaction products were run on 
an ABI 3130 Genetic Analyzer (Applied Biosystems). The 
resulting sequences were edited using CodonCodeAlign- 
er V.7 (CodonCode Corporation). Identified variable sites 
were checked by eye on the original chromatograph file in 
BioEdit (Hall 1999). Sequences were aligned in BioEdit 
(Hall 1999) with Clustal W (Thompson et al. 1994). Cyt-b 
was translated into amino acids for authentication on Ex- 



Figure I. (A) Sampling locations of Rana graeca along the species range (taken from Sunje et al. 2018). The color of each point 
indicates the haplotype inferred from the respective sample code (N = 4). The grayscale map is drawn from the NASA elevation 
model SRTMGLl (Farr & Kobrick 2000). (B) Zoomed area with sampling codes, respective locations and countries (cf Table 1). 
The baseline map is showing river flows and tributaries in light blue (www.sciencebase.gov/catalog/item/537f6c6be4b021317a87279a). 
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Table 1. Sampling codes and locations with geographic coordinates in WGS84 format. BIH - Bosnia and Herzegovina, MNE - 
Montenegro, ALB - Albania, GRC - Greece. PZC - private Zoological collection, NHMW - Natural History museum of Vienna, 
HNHM - Hungarian Natural History Museum. The code in the parenthesis, following the collection/museum abbreviation, indicates 
its registration number in the respective collection. Nb - Number of samples (total Nb = 22). H - Code of inferred haplotype. EMBL 
An - GenBank accession number. 


Code 

Location 

N 

E 

Source 

Nb 

H 

EMBL An 

CZ 

Gazin, BIH 

44.96 

15.94 

PZC: „Lelo“ (Rg9L) 

1 

1 

MH607444 

VB 

Vrelo Bosne, BIH 

43.81 

18.26 

Field work 

1 

1 

MH607446 

DARV 

Dariva, Bentbasa, BIH 

43.85 

18.45 

Field work 

1 

1 

MH607445 

ILV 

Ilovice, BIH 

43.68 

18.44 

NHMW (30428/1) 

1 

1 

MH607447 

NRT 

Neretva - Daici, BIH 

43.59 

18.04 

Field work 

1 

1 

MH607450 

RAK 

Rakitnica - Kasici, BIH 

43.55 

18.03 

Field work 

4 

1 

MH607451-54 

TUS 

TusUacka rijeka - Visocica, BIH 

43.62 

18.27 

Field work 

1 

1 

MH607460 

CJZR 

Crno jezero - Zelengora, BIH 

43.38 

18.58 

Field work 

1 

1 

MH607443 

SUH 

Suha - Cemerno, BIH 

43.23 

18.61 

NHMW (30428/5) 

1 

1 

MH607459 

SUT 

Sutjeska, BIH 

43.31 

18.66 

Field work 

3 

1 

MH607455-57 

NIK 

Niksic, MNE 

42.81 

18.88 

EMBL (Yuan et al. 2016) 

1 

1 

KX269345.1 

HRCNV 

Herceg Novi, MNE 

42.47 

18.54 

Field work 

1 

2 

MH607461 

PTR 

Petrovac, MNE 

42.20 

18.95 

Field work 

1 

3 

MH607463 

BUS 

Bushtrice, ALB 

41.90 

20.39 

HNHM (2358) 

1 

1 

MH607458 

MUR 

Murres, ALB 

41.60 

20.38 

HNHM (2357) 

1 

1 

MH607449 

LIQ 

Liqeni, ALB 

41.79 

20.19 

HNHM (2355) 

1 

1 

MH607448 

MTR 

Meteora, GRC 

39.84 

21.91 

NHMW (29072/1) 

1 

4 

MH607462 


pasy server (https://web.expasy.org/translate/). Sequences 
are deposited in the EMBL database (Table i). 

The number of haplotypes, haplotype diversity, number 
of polymorphic sites (singletons and parismony informa¬ 
tive) as well as nucleotide diversity were computed with 
DnaSp (Librado & Rozas 2009). The same software was 
used to assess Strobeck’s S, Fu and Li’s statistics and Taji- 
ma D test. Due to the fact that the number of samples from 
same sampling locations was overall too low (Table 1), it was 
not considered meaningful to define distinct populations 
which would allow additional population genetic tests. 

For the phylogenetic analyses we downloaded the single 
available R. graeca sequence longer than 500 bp available in 
EMBL Genbank database (accession number: KX269345.1, 
Yuan et al. 2016). We included the sequence of R. dalma- 
tina from this study (accession number: MH607464) in 
the analyses and specified it as outgroup following the re¬ 
sults of previous molecular analyses that included different 
Rana species (Vences et al. 2013, Yuan et al. 2016). Phylo¬ 
genetic relationships were estimated using Maximum 
Likelihood (ML) and Bayesian Inference (BI). Partition of 
the dataset was carried out in DAMBE (XiA 2001). JModel- 
Test 2.1.1 (Darriba et al. 2012) was used to infer the best 
evolutionary model under the AIC criterion for each par¬ 
tition separately and for the original alignment. To test if 
the partitioned dataset is of good quality for phylogenetic 
analysis we performed the substitution saturation test in 
DAMBE (XiA et al. 2003) for each partition separately. ML 
analyses were done in MEGA 7 (Kumar et al. 2015) using 
the bootstrap method (1,000 times) under the Tamura- 
Nei model, as suggested from the results of JModelTest. BI 


was done in MRBAYES 3.2.6 (Huelsenbeck & Ronquist 
2001); for each partition separately, we specified respective 
model parameters based on the results obtained in JModel- 
test. Bayesian analysis were done for two independent runs 
of 2 million generations, each with four Markov Ghains, 
discarding a 10% as burn-in, and sampled every 1,000 gen¬ 
erations. The tree was visualized and edited with FigTree v. 
1.3.1 (Rambaut 2009). 

Divergence and geographical distribution of the haplo¬ 
types were analyzed by means of TGS networks (Glement 
et al. 2000), using POPART 1.7 (Leight & Brian 2015). 

The amplified cyt-b was 554 bp in length. Due to miss¬ 
ing data and gaps (N = 66), 488 sites were included in the 
analysis. In total, eight sites revealed to be polymorphic; 
all were singletons with no parsimony informative sites. 
The nucleotide diversity (Pi) was 0.00149 and the average 
number of nucleotide differences (k) was 0.727. Four (4) 
distinct haplotypes were revealed in this study (Table 1, Fig. 
1 and Fig. 2) showing a haplotype diversity of 0.26 (sd = 
0.120). Strobeck’s S statistic (0.861, P = 0.207) confirmed 
that the number of haplotypes is lower or equal to four. Fu 
and Li’s statistic proved that all mutations are not selective¬ 
ly neutral (D = -3.32, P < 0.02; F = -3.19, P < 0.02), where¬ 
as the Tajima’s D suggested a recent population expansion 
(D = -2.19, P < 0.01). 

The substitution saturation test showed that the sequenc¬ 
es have experienced little substitution saturation (P = 0.000 
for all partitions), and therefore are suitable for phylo¬ 
genetic analysis. Both ML and BI phylogenetic trees of the 
haplotypes were consistent in topology and confirmed that 
R. graeca is monotypic with respect to R. dalmatina. Based 
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on the phylogenetic tree (Fig. 2), all sequences seem to be 
part of one clade but the haplotype network (Fig. 2) sug¬ 
gests that the sample from Greece might belong to another 
lineage differentiated from the rest of the populations. 

Overall, a low variation of mtDNA has been revealed in 
this study. The largest number of individuals represented a 
single haplotype that is present across the entire northern 
part of the species’ range (from BIH to Albania - Table 1, 
Figs 1-2). The sample from the southern part of the spe¬ 
cies’ range is represented by a more diverged mitochondri¬ 
al haplotype (Fig. 2) suggesting a different colonization ori¬ 
gin of these populations possibly located in Greece. Glacial 
refugia and colonization events that started from Greece 
are reported for several other taxa, such as plants (Bros 
2010), amphibians (Pabijan et al. 2014), reptiles (Lenk et 
al. 1999, PouLAKAKis et al. 2005) and birds (Griswold & 
Baker 2002, Brito 2005). The area of Montenegro harbors 
two haplotypes and is characterized by higher genetic vari¬ 
ation, which suggests that this region could have been a di¬ 
versification center for the northern populations. Ursen- 
BACHER et al. (2008) also report the area of Montenegro 
as an important diversification center explaining the phy- 
logeography of the nose horned Viper (Vipera ammodytes 
Linnaeus, 1758). The existence of multiple glacial refugia 
in the Balkans has been reported by many authors (e.g. see 
afore mentioned) but well showed in Vences et al. (2013) 
where the authors assessed these using environmental 
niche modelling. 

During the last glacial maxima the Balkans were mostly 
free of ice (Fig. iB), but the glaciation has made a marked 



Figure 2. Phylogenetic tree and haplotype network of 22 indi¬ 
viduals sampled along the species’ range. In the phylogenetic 
tree, samples with identical sequences are separated with an as¬ 
terisk (*); labels are as in Table 1 and Figure IB. Values on nodes 
correspond to number of substitution per site (the absence of a 
number in front of a sequence indicates no substitution - 0 . 000 ). 
In the haplotype network, each identified haplotype (N = 4) is a 
separate chart pie; colors correspond to countries of origin (red = 
Bosnia and Flerzegovina, green = Montenegro, purple = Albania, 
yellow = Greece); the minimum number of substitutions between 
haplotypes corresponds to the number of bars on the branches of 
the haplotype network. 


effect on the Mediterranean climate that might have been 
responsible for a retreat of brown frogs within areas of their 
distribution range (Veith et al. 2003, Sunje et al. 2018). 
Possibly, at that time, the dispersal events of R. graeca were 
ongoing to a certain extent during warm phases but were 
probably restrained by climatic oscillations and the ex¬ 
istence of several geographic barriers (e.g. mountain gla¬ 
ciers). These might have influenced a spread of populations 
from North to South and vice versa. 

A lower diversity characterized by homogenous and in¬ 
variable sequences was detected also in R. dalmatina when 
compared to R. temporaria, which is also explained by the 
more narrow distribution of the former when compared to 
the latter (Veith et al. 2003). The distribution of R. graeca 
is even narrower when compared to R. dalmatina, there¬ 
fore the observed low variation revealed in this study may 
not be surprising. A more shallow and derived lineages 
and low variation is characteristic to areas of recent range 
expansion (Veith et al. 2013) which seems also to be the 
case of R. graeca as showed by the Tajima test. Addition¬ 
ally, only limited genetic variation and a small number of 
haplotypes occur in recently (re)colonized areas due to 
range expansions that are associated with an overall re¬ 
duction of variation (Excoffier et al. 2009). Based on 
these findings, our result suggest that R. graeca is expand¬ 
ing northwards since a single haplotype is present across 
the entire northern part of the species’ range. Fast postgla¬ 
cial range expansions associated with genetic uniformity 
have been reported also for other amphibian species (e.g. 
see Kuchta & Tan 2005, Makowsky et al. 2009). The rich 
river net along the range consisting of interconnected riv¬ 
ers, streams, and tributaries (Fig. iB) seem to be facilitating 
dispersal events and expansion. Although dispersal abili¬ 
ties of R. graeca are unknown, observations far from wa¬ 
ter sources have been reported (see Bolkay 1919, Dzukic 
1968) which pose the question of the species’ plasticity. The 
registered presence of the species in unconventional habi¬ 
tats in the northern part of the range (Sunje et al. 2018) 
suggests a still ongoing expansion and the possible exist¬ 
ence of traits that favor species persistence and increase its 
potential to occupy novel habitats. By unconventional hab¬ 
itats we refer to karstic fields with (sinking) rivers, the beds 
of which are armored with grass rather than rocks, which 
is atypical when compared to rocky canyons of fast flowing 
rivers that are common to R. graeca (Dzukic 1968). The 
(adaptive) variability in ecological, morphological, behav¬ 
ioral, and physiological traits remains still to be studied to 
confirm this speculation. A wider sampling, together with 
additional genetic markers are needed to further clarify the 
evolutionary relationships within R. graeca. 
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